setwd("C:/Users/forbing36/Desktop/meta分析/Analysis")
library("meta")
library("metafor")
# install.packages("numDeriv") 
# install.packages("lme4")
# install.packages("BiasedUrn")
library(grid)
data0 <- read.csv("./matrix-OR.csv")

# TYPE：2−year Local Control ####

#修改保存路径 2y-LC-type_add
setwd("C:/Users/forbing36/Desktop/meta分析/Analysis/result2-bin/2y-LC-type_add")

## 1. meta-type####
data1<-data0[c('Author','Year',
               'events.2L.UC', #修改
               'total.2L.UC',  #修改
               'events.2L.C',  #修改
               'total.2L.C',
               'Type.UC','Type.C')]  #修改，type亚组
data1
data2<-na.omit(data1)

# which(data2$Author=="Sidiqi")
# which(data2$Author=="Franceschini")
data2<-data2[-1,]#去掉type.c和type.uc亚组不一致的行
data2<-data2[-1,]
data.Type

data.method<-data2
data00<-data2

meta0 <- metabin(events.2L.UC, #修改
                 total.2L.UC,  #修改
                 events.2L.C,  #修改
                 total.2L.C,   #修改
                 data = data.method,
                 method = "MH",
                 method.tau="DL",
                 sm="OR",
                 subgroup = Type.UC,
                 fixed= TRUE,random = TRUE,
                 studlab =paste(data00$Author,data00$Year,sep = "-"),
                 allstudies = TRUE)#,incr=0.5,allincr=TRUE,addincr=TRUE

meta1 <- metabin(events.2L.UC, #修改
                 total.2L.UC,  #修改
                 events.2L.C,  #修改
                 total.2L.C,   #修改
                 data = data.method,
                 method = "MH",
                 method.tau="DL",
                 sm="OR",
                 fixed= TRUE,random = TRUE,
                 studlab =paste(data00$Author,data00$Year,sep = "-"),
                 allstudies = TRUE)#,incr=0.5,allincr=TRUE,addincr=TRUE
summary(meta0)
output <- capture.output(summary(meta0))# 调用函数并捕获输出
cat(output, file = "meta-type.txt", sep = "\n")

# ~Method  ？？
# meta2<-metareg(meta0,~Method,intercept =TRUE)
# meta2
# meta.method <- metabin(events.2L.UC, #修改
#                        total.2L.UC,  #修改
#                        events.2L.C,  #修改
#                        total.2L.C,   #修改
#                        data = data.method,
#                        method = "MH",
#                        method.tau="DL",
#                        sm="OR",
#                        fixed= TRUE,
#                        random = TRUE,
#                        subgroup = Method,
#                        studlab =paste(data00$Author,data00$Year,sep = "-"))#,incr=0.5,allincr=TRUE,addincr=TRUE
# meta.method
# summary(meta.method)


# forest(meta1)
pdf(file="2y-LC-type.pdf", width = 13, height = 6)
forest(meta0,family="sans",fontsize=9.5,label.e="ultracentral",label.c="central",
       lwd=2,col.diamond.fixed="lightslategray",col.diamond.lines.fixed="lightslategray",
       col.diamond.random="maroon",col.diamond.lines.random="maroon",col.square="skyblue",col.study="lightslategray",
       lty.fixed=4,plotwidth="8cm",colgap.forest.left="1cm",colgap.forest.right="1cm",just.forest="right",colgap.left="0.5cm",
       colgap.right="0.5cm")
grid.text("2−year Local Control", #修改
          .53, .93, gp=gpar(cex=1))#add title
dev.off()


## 2. 漏洞图####
#漏斗图：发表偏倚
pdf(file="funel.pdf", width = 8, height = 6)
funnel(meta1,studlab = FALSE)
dev.off()

## 3. 线性图####
### peters####
metabias(meta1,method="peters",k.min = 7)#线性图：发表偏倚
output <- capture.output(metabias(meta1,method="peters",k.min = 5))
cat(output, file = "peters.txt", sep = "\n")

### egger####
metabias(meta1,method="egger",k.min=7)#线性图：发表偏倚,可不用
output <- capture.output(metabias(meta1,method="egger",k.min=7))# 调用函数并捕获输出
cat(output, file = "egger.txt", sep = "\n")

### Thompson####
metabias(meta1,method="Thompson",k.min=7)#线性图：发表偏倚
output <- capture.output(metabias(meta1,method="Thompson",k.min=7))# 调用函数并捕获输出
cat(output, file = "Thompson.txt", sep = "\n")

### linreg####
pdf(file="linreg.pdf", width = 8, height = 6)
metabias(meta1,method="linreg",plotit = TRUE,k.min=7)#线性图：发表偏倚
dev.off()
output <- capture.output(metabias(meta1,method="linreg",plotit = TRUE,k.min=7))# 调用函数并捕获输出
cat(output, file = "linreg.txt", sep = "\n")


##4. sensitive####
metainf(meta1,pooled = "fixed")
pdf(file="sensitive.pdf", width = 13, height = 6)
# forest(metainf(meta1,pooled = "fixed"),fixed=TRUE,family="sans",fontsize=10,lwd=2,col.diamond.fixed="lightslategray",col.diamond.lines.fixed="lightslategray",
#        col.diamond.random="maroon",col.diamond.lines.random="maroon",col.square="skyblue",col.study="lightslategray",
#        lty.fixed=4,plotwidth="8cm",colgap.forest.left="1cm",colgap.forest.right="1cm",just.forest="right",colgap.left="0.5cm",
#        colgap.right="0.5cm")#xlim = c(0,1),X轴范围
forest(metainf(meta1,pooled = "random"),random=TRUE,family="sans",fontsize=10,lwd=2,col.diamond.fixed="lightslategray",col.diamond.lines.fixed="lightslategray",
       col.diamond.random="maroon",col.diamond.lines.random="maroon",col.square="skyblue",col.study="lightslategray",
       lty.fixed=4,plotwidth="8cm",colgap.forest.left="1cm",colgap.forest.right="1cm",just.forest="right",colgap.left="0.5cm",
       colgap.right="0.5cm")#xlim = c(0,1),X轴范围
dev.off()

##5. trimfill####
tf1 <- trimfill(meta1,common= TRUE,random = TRUE)#com.random= FALSE
summary(tf1)
output <- capture.output(summary(tf1))# 调用函数并捕获输出
cat(output, file = "tf.txt", sep = "\n")

### tf####
pdf(file="tf.pdf", width = 8, height = 6)
# funnel(tf1)
funnel(tf1,pch=ifelse(tf1$trimfill,1,10),level = 0.95,comb.fixed= TRUE)
dev.off()

### tf-forest####
pdf(file="tf-forest.pdf", width = 13, height = 6)
forest(tf1)
dev.off()


# TYPE：Grade5 Toxicity ####

#修改保存路径 Gd5-type_add
setwd("C:/Users/forbing36/Desktop/meta分析/Analysis/result2-bin/Gd5-type_add")

## 1. meta-type####
data1<-data0[c('Author','Year',
               'events.Gd5.UC', #修改
               'total.Gd5.UC',  #修改
               'events.Gd5.C',  #修改
               'total.Gd5.C',
               'Type.UC','Type.C')]  #修改，type亚组
data1
data2<-na.omit(data1)

which(data2$Author=="Daly")
which(data2$Author=="Korzets ceder")
which(data2$Author=="Chang")
data2<-data2[-c(3,5,6),]
# data.Type

data.method<-data2
data00<-data2
meta0 <- metabin(events.Gd5.UC, #修改
                 total.Gd5.UC,  #修改
                 events.Gd5.C,  #修改
                 total.Gd5.C,   #修改
                 data = data.method,
                 method = "MH",
                 method.tau="DL",
                 sm="OR",
                 subgroup = Type.UC,
                 fixed= TRUE,random = TRUE,
                 studlab =paste(data00$Author,data00$Year,sep = "-"),
                 allstudies = TRUE)#,incr=0.5,allincr=TRUE,addincr=TRUE
summary(meta0)
output <- capture.output(summary(meta0))# 调用函数并捕获输出
cat(output, file = "meta-type.txt", sep = "\n")

meta1 <- metabin(events.Gd5.UC, #修改
                 total.Gd5.UC,  #修改
                 events.Gd5.C,  #修改
                 total.Gd5.C,   #修改
                 data = data.method,
                 method = "MH",
                 method.tau="DL",
                 sm="OR",
                 fixed= TRUE,random = TRUE,
                 studlab =paste(data00$Author,data00$Year,sep = "-"),
                 allstudies = TRUE)#,incr=0.5,allincr=TRUE,addincr=TRUE

# ~Method  ？？
# meta2<-metareg(meta0,~Method,intercept =TRUE)
# meta2
# meta.method <- metabin(events.Gd5.UC, #修改
#                        total.Gd5.UC,  #修改
#                        events.Gd5.C,  #修改
#                        total.Gd5.C,   #修改
#                        data = data.method,
#                        method = "MH",
#                        method.tau="DL",
#                        sm="OR",
#                        fixed= TRUE,
#                        random = TRUE,
#                        subgroup = Method,
#                        studlab =paste(data00$Author,data00$Year,sep = "-"))#,incr=0.5,allincr=TRUE,addincr=TRUE
# meta.method
# summary(meta.method)


# forest(meta1)
pdf(file="Gd5-type.pdf", width = 13, height = 8)
forest(meta0,family="sans",fontsize=9.5,label.e="ultracentral",label.c="central",
       lwd=2,col.diamond.fixed="lightslategray",col.diamond.lines.fixed="lightslategray",
       col.diamond.random="maroon",col.diamond.lines.random="maroon",col.square="skyblue",col.study="lightslategray",
       lty.fixed=4,plotwidth="8cm",colgap.forest.left="1cm",colgap.forest.right="1cm",just.forest="right",colgap.left="0.5cm",
       colgap.right="0.5cm")
grid.text("Grade 5 Toxicity", #修改
          .53, .92, gp=gpar(cex=1))#add title
dev.off()


## 2. 漏洞图####
#漏斗图：发表偏倚
pdf(file="funel.pdf", width = 8, height = 6)
funnel(meta1,studlab = FALSE)
dev.off()

## 3. 线性图####
### peters####
metabias(meta1,method="peters",k.min = 9)#线性图：发表偏倚
output <- capture.output(metabias(meta1,method="peters",k.min = 9))
cat(output, file = "peters.txt", sep = "\n")

### egger####
metabias(meta1,method="egger",k.min=9)#线性图：发表偏倚,可不用
output <- capture.output(metabias(meta1,method="egger",k.min=7))# 调用函数并捕获输出
cat(output, file = "egger.txt", sep = "\n")

### Thompson####
metabias(meta1,method="Thompson",k.min=9)#线性图：发表偏倚
output <- capture.output(metabias(meta1,method="Thompson",k.min=7))# 调用函数并捕获输出
cat(output, file = "Thompson.txt", sep = "\n")

### linreg####
pdf(file="linreg.pdf", width = 8, height = 6)
metabias(meta1,method="linreg",plotit = TRUE,k.min=9)#线性图：发表偏倚
dev.off()
output <- capture.output(metabias(meta1,method="linreg",plotit = TRUE,k.min=9))# 调用函数并捕获输出
cat(output, file = "linreg.txt", sep = "\n")


##4. sensitive####
metainf(meta1,pooled = "fixed")
pdf(file="sensitive.pdf", width = 13, height = 6)
# forest(metainf(meta1,pooled = "fixed"),fixed=TRUE,family="sans",fontsize=10,lwd=2,col.diamond.fixed="lightslategray",col.diamond.lines.fixed="lightslategray",
#        col.diamond.random="maroon",col.diamond.lines.random="maroon",col.square="skyblue",col.study="lightslategray",
#        lty.fixed=4,plotwidth="8cm",colgap.forest.left="1cm",colgap.forest.right="1cm",just.forest="right",colgap.left="0.5cm",
#        colgap.right="0.5cm")#xlim = c(0,1),X轴范围
forest(metainf(meta1,pooled = "random"),random=TRUE,family="sans",fontsize=10,lwd=2,col.diamond.fixed="lightslategray",col.diamond.lines.fixed="lightslategray",
       col.diamond.random="maroon",col.diamond.lines.random="maroon",col.square="skyblue",col.study="lightslategray",
       lty.fixed=4,plotwidth="8cm",colgap.forest.left="1cm",colgap.forest.right="1cm",just.forest="right",colgap.left="0.5cm",
       colgap.right="0.5cm")#xlim = c(0,1),X轴范围
dev.off()

##5. trimfill####
tf1 <- trimfill(meta1,common= TRUE,random = TRUE)#com.random= FALSE
summary(tf1)
output <- capture.output(summary(tf1))# 调用函数并捕获输出
cat(output, file = "tf.txt", sep = "\n")

### tf####
pdf(file="tf.pdf", width = 8, height = 6)
# funnel(tf1)
funnel(tf1,pch=ifelse(tf1$trimfill,1,10),level = 0.95,comb.fixed= TRUE)
dev.off()

### tf-forest####
pdf(file="tf-forest.pdf", width = 13, height = 6)
forest(tf1)
dev.off()

##6. sparse data analysis####
# data1<-data0[c('Author',
#                'Year',
#                'events.Gd5.UC', #修改
#                'total.Gd5.UC',  #修改
#                'events.Gd5.C',  #修改
#                'total.Gd5.C')]  #修改
# data1<-data0[c('Author','Year','Type.UC','Type.C','events.Gd3.UC','total.Gd3.UC','events.Gd3.C','total.Gd3.C')]
# data1
# data2<-na.omit(data1)
# data2
# data.BED<-data2
# which(data2$Author=="Regnery")
# which(data2$Author=="Franceschini")
# data2<-data2[-12,]
# data2<-data2[-8,]
# data.Type
# data00<-data.BED
data00
meta1<-rma.glmm(measure = "OR",
                ai=events.Gd5.UC, #修改
                n1i=total.Gd5.UC, #修改
                ci=events.Gd5.C,  #修改
                n2i=total.Gd5.C,  #修改
                data = data00,
                slab =paste(data00$Author,data00$Year,sep = "-"),
                drop00=FALSE,method="FE",model="CM.EL")
###meta-GLMM####
print(meta1,digits=3) #GLMM方法
output <- capture.output(print(meta1,digits=3) )# 调用函数并捕获输出
cat(output, file = "meta-GLMM.txt", sep = "\n")

###meta-GLMM-predict####
predict(meta1,transf=exp,digits=3)
summary(meta1)
output <- capture.output(summary(meta1))# 调用函数并捕获输出
cat(output, file = "meta-GLMM-predict.txt", sep = "\n")

### meta-GLMM-forest####
# forest(meta1)
pdf(file="meta-GLMN.pdf", width = 13, height = 6)
forest(meta1,xlim = c(-20, 10),at = log(c(0.05, 0.25, 1, 4,8)),atransf = exp,
       ilab = cbind(data00$events.Gd5.UC, data00$total.Gd5.UC, data00$events.Gd5.C, data00$total.Gd5.C),ilab.xpos = c(-13, -11, -9, -7),
       cex=1, header="Author-Year",mlab="",psize = 1)
op <- par(cex=1, font=2)
text(c(-13, -11, -9, -7), 16.5, c("Events", "Total", "Events", "Total"),cex =1)
text(c(-12,-8), 17, c("Ultracentral", "Central"),cex = 1)
text(0,17,"Grade5 Toxicity",cex = 1)#修改
par(op)

#添加 Q 值、自由度、p 值等统计量
text(-20, -1, pos = 4, cex = 0.8, bquote(paste("FE Model (Q.WLD(df = 14)=",
                                               .(formatC(meta1$QE.Wld,digits=2,format="f")),
                                               ",p=",.(formatC(meta1$QEp.Wld,digits=2,format="f")),
                                               ";Q.LRT(df = 14)=",.(formatC(meta1$QE.LRT,digits=2,format="f")),
                                               ",p=",.(formatC(meta1$QEp.LRT,digits=2,format="f")),";",I^2,"=",
                                               .(formatC(meta1$I2,digits=1,format="f")),"%)")))
dev.off()

#？？？？
# meta1<-metabin(event.e = events.Gd5.UC,
#                n.e = total.Gd5.UC,
#                event.c = events.Gd5.C, 
#                n.c = total.Gd5.C,
#                studlab= paste(data00$Author,data00$Year,sep = "-"), data =data00,sm="OR",
#                method = "MH", method.tau= "DL",allstudies=TRUE)

##7. normal data __未使用####
data1<-data0[c('Author','Year','events.Gd3.UC','total.Gd3.UC','events.Gd3.C','total.Gd3.C')]
data1
data2<-na.omit(data1)
data2
data.PTV<-data2
which(data2$Author=="Sidiqi")
which(data2$Author=="Franceschini")
which(data2$Author=="Regnery")
data2<-data2[-6,]
data2<-data2[-1,]
data.PTV
data.PTV<-data.PTV[order(data.PTV$PTV.UC),]
data00<-data.PTV
data00
meta0 <- metabin(events.Gd3.UC,total.Gd3.UC,events.Gd3.C,total.Gd3.C,data = data00,method = "MH",method.tau="DL",sm="OR",fixed= TRUE,random = TRUE,studlab =paste(data.PTV$Author,data.PTV$Year,sep = "-"),allstudies = TRUE)#,incr=0.5,allincr=TRUE,addincr=TRUE
summary(meta0)
meta2<-metareg(meta0,~PTV.UC,intercept =TRUE)
meta2


meta1 <-metabin(events.1L.UC,total.1L.UC,events.1L.C,total.1L.C,data = data.PTV,method = "MH",method.tau="DL",sm="OR",fixed= TRUE,random = TRUE,studlab =paste(data.PTV$Author,data.PTV$Year,sep = "-"))
meta1
summary(meta1)
forest(meta1)
forest(meta0,family="sans",fontsize=10,label.e="ultracentral",label.c="central",
       lwd=2,col.diamond.fixed="lightslategray",col.diamond.lines.fixed="lightslategray",
       col.diamond.random="maroon",col.diamond.lines.random="maroon",col.square="skyblue",col.study="lightslategray",
       lty.fixed=4,plotwidth="8cm",colgap.forest.left="1cm",colgap.forest.right="1cm",just.forest="right",colgap.left="0.5cm",
       colgap.right="0.5cm")
grid.text("2-year Local Control", .55, .73, gp=gpar(cex=1))#add title
grid.text("Grade3-5 Toxicity", .53, .76, gp=gpar(cex=1))#add title
grid.text("5-year Overall Survival", .55, .76, gp=gpar(cex=1))#add title
dev.off()
funnel(meta0,studlab = FALSE,cex=1.5)
metabias(meta0,method="peters")
metabias(meta0,method="egger")
metabias(meta0,method="Thompson",k.min=10)
metabias(meta0,method="linreg",plotit = TRUE,k.min=7)
jpeg("Sensitivity analysis results.jpeg",height = 1300,width = 800)
metainf(meta0,pooled = "fixed")
forest(metainf(meta0,pooled = "fixed"),fixed =TRUE,family="sans",fontsize=10,lwd=2,col.diamond.fixed="lightslategray",col.diamond.lines.fixed="lightslategray",
       col.diamond.random="maroon",col.diamond.lines.random="maroon",col.square="skyblue",col.study="lightslategray",
       lty.fixed=4,plotwidth="8cm",colgap.forest.left="1cm",colgap.forest.right="1cm",just.forest="right",colgap.left="0.5cm",
       colgap.right="0.5cm")#xlim = c(0,1),X轴范围
dev.off()
tf1 <- trimfill(meta0,common= TRUE,random = TRUE)#com.random= FALSE
summary(tf1)
funnel(tf1)
funnel(tf1,pch=ifelse(tf1$trimfill,1,10),level = 0.95,comb.fixed= TRUE,cex=1.5)
forest(tf1)
forest(tf1,family="sans",fontsize=10,label.e="ultracentral",label.c="central",test.subgroup = TRUE,
       lwd=2,col.diamond.fixed="lightslategray",col.diamond.lines.fixed="lightslategray",
       col.diamond.random="maroon",col.diamond.lines.random="maroon",col.square="skyblue",col.study="lightslategray",
       lty.fixed=4,plotwidth="8cm",colgap.forest.left="1cm",colgap.forest.right="1cm",just.forest="right",colgap.left="0.5cm",
       colgap.right="0.5cm")


#举例
meta1<-metabin(event.e = Events.UC,n.e = Total.UC,event.c = Events.C, n.c = Total.C,
               studlab= paste(data$Author,data$Year,sep = "-"), data =data,sm="OR",
               method = "GLMM", model.glmm ="CM.EL",method.tau= "ML",nAGQ=1)
#analysis of subgroup
forest(meta1,family="sans",fontsize=9.5,label.e="ultracentral",label.c="central",
       lwd=2,col.diamond.fixed="lightslategray",col.diamond.lines.fixed="lightslategray",
       col.diamond.random="maroon",col.diamond.lines.random="maroon",col.square="skyblue",col.study="lightslategray",
       lty.fixed=4,plotwidth="8cm",colgap.forest.left="1cm",colgap.forest.right="1cm",just.forest="right",colgap.left="0.5cm",
       colgap.right="0.5cm",test.subgroup=TRUE,resid.hetstat=TRUE,allstudies = FALSE)










#PTV：2−year Overall Survival####
#修改保存路径2y-OS-PTV_add
setwd("C:/Users/forbing36/Desktop/meta分析/Analysis/result2-bin/2y-OS-PTV_add")

## 1. meta-type####
data1<-data0[c('Author','Year',
               'PTV.UC',
               'PTV.C',
               'events.2O.UC', #修改
               'total.2O.UC',  #修改
               'events.2O.C',  #修改
               'total.2O.C')]  #修改
data1
data2<-na.omit(data1)
data2
data.PTV<-data2
# which(data2$Author=="Sidiqi")
# which(data2$Author=="Franceschini")
# which(data2$Author=="Regnery")
# data2<-data2[-6,]
# data2<-data2[-1,]
data.PTV
data.PTV<-data.PTV[order(data.PTV$PTV.UC),]#排序
data00<-data.PTV
data00
meta0 <- metabin(events.2O.UC, #修改
                 total.2O.UC,  #修改
                 events.2O.C,  #修改
                 total.2O.C,   #修改
                 data = data.PTV,
                 method = "MH",
                 method.tau="DL",
                 sm="OR",
                 fixed= TRUE,
                 random = TRUE,
                 studlab =paste(data.PTV$Author,data.PTV$Year,sep = "-"))#,incr=0.5,allincr=TRUE,addincr=TRUE
summary(meta0)

output <- capture.output(summary(meta0))# 调用函数并捕获输出
cat(output, file = "meta-type.txt", sep = "\n")

meta2<-metareg(meta0,~PTV.UC,intercept =TRUE)
meta2
meta.PTV<- metabin(events.2O.UC,total.2O.UC,events.2O.C,total.2O.C,data = data.PTV,method = "MH",method.tau="DL",sm="OR",fixed= TRUE,random = TRUE,
                   subgroup = ifelse(PTV.UC<62.5,"PTV.UC<62.5","PTV.UC>=62.5"),studlab =paste(data.PTV$Author,data.PTV$Year,sep = "-"))#,incr=0.5,allincr=TRUE,addincr=TRUE
meta.PTV
summary(meta.PTV)

output <- capture.output(summary(meta.PTV))# 调用函数并捕获输出
cat(output, file = "meta_PTV.txt", sep = "\n")

# forest(meta0)
pdf(file="2y-OS-PTV.pdf", width = 13, height = 6)
forest(meta.PTV,family="sans",fontsize=10,label.e="ultracentral",label.c="central",test.subgroup = TRUE, 
       lwd=2,col.diamond.fixed="lightslategray",col.diamond.lines.fixed="lightslategray",
       col.diamond.random="maroon",col.diamond.lines.random="maroon",col.square="skyblue",col.study="lightslategray",
       lty.fixed=4,plotwidth="8cm",colgap.forest.left="1cm",colgap.forest.right="1cm",just.forest="right",colgap.left="0.5cm",
       colgap.right="0.5cm",print.subgroup.name = FALSE)
grid.text("2-year Overall Survival", .52, .96, gp=gpar(cex=1))#add title
# grid.text("2-year Local Control", .52, .86, gp=gpar(cex=1))
dev.off()

##2. funnel####
pdf(file="funnel.pdf", width = 8, height = 6)
funnel(meta.PTV,studlab = FALSE,cex=1.5)
dev.off()

# metabias(meta1,method="peters")
# metabias(meta1,method="egger",k.min=7)
# metabias(meta1,method="Thompson",k.min=10)
# metabias(meta1,method="linreg",plotit = TRUE,k.min=7)
# jpeg("Sensitivity analysis results.jpeg",height = 1300,width = 800)
# metainf(meta0,pooled = "fixed")


##3. forest####
pdf(file="forest.pdf", width = 13, height = 6)
forest(metainf(meta0,pooled = "fixed"),fixed=TRUE,family="sans",fontsize=9.5,lwd=2,col.diamond.fixed="lightslategray",col.diamond.lines.fixed="lightslategray",
       col.diamond.random="maroon",col.diamond.lines.random="maroon",col.square="skyblue",col.study="lightslategray",
       lty.fixed=4,plotwidth="8cm",colgap.forest.left="1cm",colgap.forest.right="1cm",just.forest="right",colgap.left="0.5cm",
       colgap.right="0.5cm")#xlim = c(0,1),X轴范围
dev.off()


##4. tf####
tf1 <- trimfill(meta0,common= TRUE,random = TRUE)#com.random= FALSE
summary(tf1)
output <- capture.output(summary(tf1))# 调用函数并捕获输出
cat(output, file = "tf.txt", sep = "\n")
pdf(file="tf.pdf", width = 8, height = 6)
# funnel(tf1)
funnel(tf1,pch=ifelse(tf1$trimfill,1,10),level = 0.95,comb.fixed= TRUE,cex=1.5)
dev.off()

##5. tf-forest####
pdf(file="tf-forest.pdf", width = 13, height = 6)
forest(tf1,family="sans",fontsize=10,label.e="ultracentral",label.c="central",test.subgroup = TRUE,
       lwd=2,col.diamond.fixed="lightslategray",col.diamond.lines.fixed="lightslategray",
       col.diamond.random="maroon",col.diamond.lines.random="maroon",col.square="skyblue",col.study="lightslategray",
       lty.fixed=4,plotwidth="8cm",colgap.forest.left="1cm",colgap.forest.right="1cm",just.forest="right",colgap.left="0.5cm",
       colgap.right="0.5cm")
dev.off()


#add####
meta1 <- meta0
## 3. 线性图####
### peters####
metabias(meta1,method="peters",k.min =8 )#线性图：发表偏倚
output <- capture.output(metabias(meta1,method="peters",k.min = 8))
cat(output, file = "peters.txt", sep = "\n")

### egger####
metabias(meta1,method="egger",k.min=8)#线性图：发表偏倚,可不用
output <- capture.output(metabias(meta1,method="egger",k.min=8))# 调用函数并捕获输出
cat(output, file = "egger.txt", sep = "\n")

### Thompson####
metabias(meta1,method="Thompson",k.min=8)#线性图：发表偏倚
output <- capture.output(metabias(meta1,method="Thompson",k.min=8))# 调用函数并捕获输出
cat(output, file = "Thompson.txt", sep = "\n")

### linreg####
pdf(file="linreg.pdf", width = 8, height = 6)
metabias(meta1,method="linreg",plotit = TRUE,k.min=7)#线性图：发表偏倚
dev.off()
output <- capture.output(metabias(meta1,method="linreg",plotit = TRUE,k.min=8))# 调用函数并捕获输出
cat(output, file = "linreg.txt", sep = "\n")



#PTV：1−year Overall Survival add_20231111####
#修改保存路径2y-OS-PTV_add
setwd("C:/Users/forbing36/Desktop/meta分析/Analysis/result2-bin/1y-OS-PTV_add")

## 1. meta-type####
data1<-data0[c('Author','Year',
               'PTV.UC',
               'PTV.C',
               'events.1O.UC', #修改
               'total.1O.UC',  #修改
               'events.1O.C',  #修改
               'total.1O.C')]  #修改
data1
data2<-na.omit(data1)
data2
data.PTV<-data2
# which(data2$Author=="Sidiqi")
# which(data2$Author=="Franceschini")
# which(data2$Author=="Regnery")
# data2<-data2[-6,]
# data2<-data2[-1,]
data.PTV
data.PTV<-data.PTV[order(data.PTV$PTV.UC),]#排序
data00<-data.PTV
data00
meta0 <- metabin(events.1O.UC, #修改
                 total.1O.UC,  #修改
                 events.1O.C,  #修改
                 total.1O.C,   #修改
                 data = data.PTV,
                 method = "MH",
                 method.tau="DL",
                 sm="OR",
                 fixed= TRUE,
                 random = TRUE,
                 studlab =paste(data.PTV$Author,data.PTV$Year,sep = "-"))#,incr=0.5,allincr=TRUE,addincr=TRUE
summary(meta0)

output <- capture.output(summary(meta0))# 调用函数并捕获输出
cat(output, file = "meta-type.txt", sep = "\n")

meta2<-metareg(meta0,~PTV.UC,intercept =TRUE)
meta2
meta.PTV<- metabin(events.1O.UC,total.1O.UC,events.1O.C,total.1O.C,data = data.PTV,method = "MH",method.tau="DL",sm="OR",fixed= TRUE,random = TRUE,
                   subgroup = ifelse(PTV.UC<62.5,"PTV.UC<62.5","PTV.UC>=62.5"),studlab =paste(data.PTV$Author,data.PTV$Year,sep = "-"))#,incr=0.5,allincr=TRUE,addincr=TRUE
meta.PTV
summary(meta.PTV)

output <- capture.output(summary(meta.PTV))# 调用函数并捕获输出
cat(output, file = "meta_PTV.txt", sep = "\n")

# forest(meta0)
pdf(file="1y-OS-PTV.pdf", width = 13, height = 6)
forest(meta.PTV,family="sans",fontsize=10,label.e="ultracentral",label.c="central",test.subgroup = TRUE, 
       lwd=2,col.diamond.fixed="lightslategray",col.diamond.lines.fixed="lightslategray",
       col.diamond.random="maroon",col.diamond.lines.random="maroon",col.square="skyblue",col.study="lightslategray",
       lty.fixed=4,plotwidth="8cm",colgap.forest.left="1cm",colgap.forest.right="1cm",just.forest="right",colgap.left="0.5cm",
       colgap.right="0.5cm",print.subgroup.name = FALSE)
grid.text("1-year Overall Survival", .52, .96, gp=gpar(cex=1))#add title
# grid.text("2-year Local Control", .52, .86, gp=gpar(cex=1))
dev.off()

##2. funnel####
pdf(file="funnel.pdf", width = 8, height = 6)
funnel(meta.PTV,studlab = FALSE,cex=1.5)
dev.off()

# metabias(meta1,method="peters")
# metabias(meta1,method="egger",k.min=7)
# metabias(meta1,method="Thompson",k.min=10)
# metabias(meta1,method="linreg",plotit = TRUE,k.min=7)
# jpeg("Sensitivity analysis results.jpeg",height = 1300,width = 800)
# metainf(meta0,pooled = "fixed")


##3. forest####
pdf(file="forest.pdf", width = 13, height = 6)
forest(metainf(meta0,pooled = "fixed"),fixed=TRUE,family="sans",fontsize=9.5,lwd=2,col.diamond.fixed="lightslategray",col.diamond.lines.fixed="lightslategray",
       col.diamond.random="maroon",col.diamond.lines.random="maroon",col.square="skyblue",col.study="lightslategray",
       lty.fixed=4,plotwidth="8cm",colgap.forest.left="1cm",colgap.forest.right="1cm",just.forest="right",colgap.left="0.5cm",
       colgap.right="0.5cm")#xlim = c(0,1),X轴范围
dev.off()


##4. tf####
tf1 <- trimfill(meta0,common= TRUE,random = TRUE)#com.random= FALSE
summary(tf1)
output <- capture.output(summary(tf1))# 调用函数并捕获输出
cat(output, file = "tf.txt", sep = "\n")
pdf(file="tf.pdf", width = 8, height = 6)
# funnel(tf1)
funnel(tf1,pch=ifelse(tf1$trimfill,1,10),level = 0.95,comb.fixed= TRUE,cex=1.5)
dev.off()

##5. tf-forest####
pdf(file="tf-forest.pdf", width = 13, height = 6)
forest(tf1,family="sans",fontsize=10,label.e="ultracentral",label.c="central",test.subgroup = TRUE,
       lwd=2,col.diamond.fixed="lightslategray",col.diamond.lines.fixed="lightslategray",
       col.diamond.random="maroon",col.diamond.lines.random="maroon",col.square="skyblue",col.study="lightslategray",
       lty.fixed=4,plotwidth="8cm",colgap.forest.left="1cm",colgap.forest.right="1cm",just.forest="right",colgap.left="0.5cm",
       colgap.right="0.5cm")
dev.off()


#add####
meta1 <- meta0
## 3. 线性图####
### peters####
metabias(meta1,method="peters",k.min =8 )#线性图：发表偏倚
output <- capture.output(metabias(meta1,method="peters",k.min = 8))
cat(output, file = "peters.txt", sep = "\n")

### egger####
metabias(meta1,method="egger",k.min=8)#线性图：发表偏倚,可不用
output <- capture.output(metabias(meta1,method="egger",k.min=8))# 调用函数并捕获输出
cat(output, file = "egger.txt", sep = "\n")

### Thompson####
metabias(meta1,method="Thompson",k.min=8)#线性图：发表偏倚
output <- capture.output(metabias(meta1,method="Thompson",k.min=8))# 调用函数并捕获输出
cat(output, file = "Thompson.txt", sep = "\n")

### linreg####
pdf(file="linreg.pdf", width = 8, height = 6)
metabias(meta1,method="linreg",plotit = TRUE,k.min=7)#线性图：发表偏倚
dev.off()
output <- capture.output(metabias(meta1,method="linreg",plotit = TRUE,k.min=8))# 调用函数并捕获输出
cat(output, file = "linreg.txt", sep = "\n")




#PTV：Grade 3-5 Toxicity####
#修改保存路径Gd3-PTV_add
setwd("C:/Users/forbing36/Desktop/meta分析/Analysis/result2-bin/Gd3-PTV_add")

## 1. meta-type####
data1<-data0[c('Author','Year',
               'PTV.UC',
               'PTV.C',
               'events.Gd3.UC', #修改
               'total.Gd3.UC',  #修改
               'events.Gd3.C',  #修改
               'total.Gd3.C')]  #修改
data1
data2<-na.omit(data1)
data2
data.PTV<-data2
# which(data2$Author=="Sidiqi")
# which(data2$Author=="Franceschini")
# which(data2$Author=="Regnery")
# data2<-data2[-6,]
# data2<-data2[-1,]
data.PTV
data.PTV<-data.PTV[order(data.PTV$PTV.UC),]#排序
data00<-data.PTV
data00
meta0 <- metabin(events.Gd3.UC, #修改
                 total.Gd3.UC,  #修改
                 events.Gd3.C,  #修改
                 total.Gd3.C,   #修改
                 data = data.PTV,
                 method = "MH",
                 method.tau="DL",
                 sm="OR",
                 fixed= TRUE,
                 random = TRUE,
                 studlab =paste(data.PTV$Author,data.PTV$Year,sep = "-"))#,incr=0.5,allincr=TRUE,addincr=TRUE
summary(meta0)
output <- capture.output(summary(meta0))# 调用函数并捕获输出
cat(output, file = "meta-type.txt", sep = "\n")

meta2<-metareg(meta0,~PTV.UC,intercept =TRUE)
meta2
meta.PTV<- metabin(events.Gd3.UC,total.Gd3.UC,events.Gd3.C,total.Gd3.C,data = data.PTV,method = "MH",method.tau="DL",sm="OR",fixed= TRUE,random = TRUE,
                   subgroup = ifelse(PTV.UC<62.5,"PTV.UC<62.5","PTV.UC>=62.5"),studlab =paste(data.PTV$Author,data.PTV$Year,sep = "-"))#,incr=0.5,allincr=TRUE,addincr=TRUE
meta.PTV
summary(meta.PTV)
output <- capture.output(summary(meta.PTV))# 调用函数并捕获输出
cat(output, file = "meta_PTV.txt", sep = "\n")
# forest(meta0)
pdf(file="Gd3-PTV.pdf", width = 13, height = 6)
forest(meta.PTV,family="sans",fontsize=10,label.e="ultracentral",label.c="central",test.subgroup = TRUE, 
       lwd=2,col.diamond.fixed="lightslategray",col.diamond.lines.fixed="lightslategray",
       col.diamond.random="maroon",col.diamond.lines.random="maroon",col.square="skyblue",col.study="lightslategray",
       lty.fixed=4,plotwidth="8cm",colgap.forest.left="1cm",colgap.forest.right="1cm",just.forest="right",colgap.left="0.5cm",
       colgap.right="0.5cm",print.subgroup.name = FALSE)
grid.text("Grade 3-5 Toxicity", .52, .96, gp=gpar(cex=1))#add title
# grid.text("2-year Local Control", .52, .86, gp=gpar(cex=1))
dev.off()

##2. funnel####
pdf(file="funnel.pdf", width = 8, height = 6)
funnel(meta.PTV,studlab = FALSE,cex=1.5)
dev.off()

# metabias(meta1,method="peters")
# metabias(meta1,method="egger",k.min=7)
# metabias(meta1,method="Thompson",k.min=10)
# metabias(meta1,method="linreg",plotit = TRUE,k.min=7)
# jpeg("Sensitivity analysis results.jpeg",height = 1300,width = 800)
# metainf(meta0,pooled = "fixed")


##3. forest####
pdf(file="forest.pdf", width = 13, height = 6)
forest(metainf(meta0,pooled = "fixed"),fixed=TRUE,family="sans",fontsize=9.5,lwd=2,col.diamond.fixed="lightslategray",col.diamond.lines.fixed="lightslategray",
       col.diamond.random="maroon",col.diamond.lines.random="maroon",col.square="skyblue",col.study="lightslategray",
       lty.fixed=4,plotwidth="8cm",colgap.forest.left="1cm",colgap.forest.right="1cm",just.forest="right",colgap.left="0.5cm",
       colgap.right="0.5cm")#xlim = c(0,1),X轴范围
dev.off()


##4. tf####
tf1 <- trimfill(meta0,common= TRUE,random = TRUE)#com.random= FALSE
summary(tf1)
output <- capture.output(summary(tf1))# 调用函数并捕获输出
cat(output, file = "tf1.txt", sep = "\n")

pdf(file="tf.pdf", width = 8, height = 6)
# funnel(tf1)
funnel(tf1,pch=ifelse(tf1$trimfill,1,10),level = 0.95,comb.fixed= TRUE,cex=1.5)
dev.off()

##5. tf-forest####
pdf(file="tf-forest.pdf", width = 13, height = 6)
forest(tf1,family="sans",fontsize=10,label.e="ultracentral",label.c="central",test.subgroup = TRUE,
       lwd=2,col.diamond.fixed="lightslategray",col.diamond.lines.fixed="lightslategray",
       col.diamond.random="maroon",col.diamond.lines.random="maroon",col.square="skyblue",col.study="lightslategray",
       lty.fixed=4,plotwidth="8cm",colgap.forest.left="1cm",colgap.forest.right="1cm",just.forest="right",colgap.left="0.5cm",
       colgap.right="0.5cm")
dev.off()

#add####
meta1 <- meta0
## 3. 线性图####
### peters####
metabias(meta1,method="peters",k.min =1 )#线性图：发表偏倚
output <- capture.output(metabias(meta1,method="peters",k.min = 1))
cat(output, file = "peters.txt", sep = "\n")

### egger####
metabias(meta1,method="egger",k.min=1)#线性图：发表偏倚,可不用
output <- capture.output(metabias(meta1,method="egger",k.min=1))# 调用函数并捕获输出
cat(output, file = "egger.txt", sep = "\n")

### Thompson####
metabias(meta1,method="Thompson",k.min=1)#线性图：发表偏倚
output <- capture.output(metabias(meta1,method="Thompson",k.min=1))# 调用函数并捕获输出
cat(output, file = "Thompson.txt", sep = "\n")

### linreg####
pdf(file="linreg.pdf", width = 8, height = 6)
metabias(meta1,method="linreg",plotit = TRUE,k.min=1)#线性图：发表偏倚
dev.off()
output <- capture.output(metabias(meta1,method="linreg",plotit = TRUE,k.min=1))# 调用函数并捕获输出
cat(output, file = "linreg.txt", sep = "\n")



#PTV：Grade 5 Toxicity####
#修改保存路径Gd5-PTV_add
setwd("C:/Users/forbing36/Desktop/meta分析/Analysis/result2-bin/Gd5-PTV_add")

## 1. meta-type####
data1<-data0[c('Author','Year',
               'PTV.UC',
               'PTV.C',
               'events.Gd5.UC', #修改
               'total.Gd5.UC',  #修改
               'events.Gd5.C',  #修改
               'total.Gd5.C')]  #修改
data1
data2<-na.omit(data1)
data2
data.PTV<-data2
# which(data2$Author=="Sidiqi")
# which(data2$Author=="Franceschini")
# which(data2$Author=="Regnery")
# data2<-data2[-6,]
# data2<-data2[-1,]
data.PTV
data.PTV<-data.PTV[order(data.PTV$PTV.UC),]#排序
data00<-data.PTV
data00
meta0 <- metabin(events.Gd5.UC, #修改
                 total.Gd5.UC,  #修改
                 events.Gd5.C,  #修改
                 total.Gd5.C,   #修改
                 data = data.PTV,
                 method = "MH",
                 method.tau="DL",
                 sm="OR",
                 fixed= TRUE,
                 random = TRUE,
                 studlab =paste(data.PTV$Author,data.PTV$Year,sep = "-"))#,incr=0.5,allincr=TRUE,addincr=TRUE
summary(meta0)
output <- capture.output(summary(meta0))# 调用函数并捕获输出
cat(output, file = "meta-type.txt", sep = "\n")


meta2<-metareg(meta0,~PTV.UC,intercept =TRUE)
meta2
meta.PTV<- metabin(events.Gd5.UC,total.Gd5.UC,events.Gd5.C,total.Gd5.C,data = data.PTV,method = "MH",method.tau="DL",sm="OR",fixed= TRUE,random = TRUE,
                   subgroup = ifelse(PTV.UC<62.5,"PTV.UC<62.5","PTV.UC>=62.5"),studlab =paste(data.PTV$Author,data.PTV$Year,sep = "-"))#,incr=0.5,allincr=TRUE,addincr=TRUE
meta.PTV
summary(meta.PTV)
output <- capture.output(summary(meta.PTV))# 调用函数并捕获输出
cat(output, file = "meta-PTV.txt", sep = "\n")

# forest(meta0)
pdf(file="Gd5-PTV.pdf", width = 13, height = 6)
forest(meta.PTV,family="sans",fontsize=10,label.e="ultracentral",label.c="central",test.subgroup = TRUE, 
       lwd=2,col.diamond.fixed="lightslategray",col.diamond.lines.fixed="lightslategray",
       col.diamond.random="maroon",col.diamond.lines.random="maroon",col.square="skyblue",col.study="lightslategray",
       lty.fixed=4,plotwidth="8cm",colgap.forest.left="1cm",colgap.forest.right="1cm",just.forest="right",colgap.left="0.5cm",
       colgap.right="0.5cm",print.subgroup.name = FALSE)
grid.text("Grade 5 Toxicity", .52, .96, gp=gpar(cex=1))#add title
# grid.text("2-year Local Control", .52, .86, gp=gpar(cex=1))
dev.off()

##2. funnel####
pdf(file="funnel.pdf", width = 8, height = 6)
funnel(meta.PTV,studlab = FALSE,cex=1.5)
dev.off()

# metabias(meta1,method="peters")
# metabias(meta1,method="egger",k.min=7)
# metabias(meta1,method="Thompson",k.min=10)
# metabias(meta1,method="linreg",plotit = TRUE,k.min=7)
# jpeg("Sensitivity analysis results.jpeg",height = 1300,width = 800)
# metainf(meta0,pooled = "fixed")


##3. forest####
pdf(file="forest.pdf", width = 13, height = 6)
forest(metainf(meta0,pooled = "fixed"),fixed=TRUE,family="sans",fontsize=9.5,lwd=2,col.diamond.fixed="lightslategray",col.diamond.lines.fixed="lightslategray",
       col.diamond.random="maroon",col.diamond.lines.random="maroon",col.square="skyblue",col.study="lightslategray",
       lty.fixed=4,plotwidth="8cm",colgap.forest.left="1cm",colgap.forest.right="1cm",just.forest="right",colgap.left="0.5cm",
       colgap.right="0.5cm")#xlim = c(0,1),X轴范围
dev.off()


##4. tf####
tf1 <- trimfill(meta0,common= TRUE,random = TRUE)#com.random= FALSE
summary(tf1)
output <- capture.output(summary(tf1))# 调用函数并捕获输出
cat(output, file = "tf1.txt", sep = "\n")
pdf(file="tf.pdf", width = 8, height = 6)
# funnel(tf1)
funnel(tf1,pch=ifelse(tf1$trimfill,1,10),level = 0.95,comb.fixed= TRUE,cex=1.5)
dev.off()

##5. tf-forest####
pdf(file="tf-forest.pdf", width = 13, height = 6)
forest(tf1,family="sans",fontsize=10,label.e="ultracentral",label.c="central",test.subgroup = TRUE,
       lwd=2,col.diamond.fixed="lightslategray",col.diamond.lines.fixed="lightslategray",
       col.diamond.random="maroon",col.diamond.lines.random="maroon",col.square="skyblue",col.study="lightslategray",
       lty.fixed=4,plotwidth="8cm",colgap.forest.left="1cm",colgap.forest.right="1cm",just.forest="right",colgap.left="0.5cm",
       colgap.right="0.5cm")
dev.off()

#add####
meta1 <- meta0
## 3. 线性图####
### peters####
metabias(meta1,method="peters",k.min =1 )#线性图：发表偏倚
output <- capture.output(metabias(meta1,method="peters",k.min = 1))
cat(output, file = "peters.txt", sep = "\n")

### egger####
metabias(meta1,method="egger",k.min=1)#线性图：发表偏倚,可不用
output <- capture.output(metabias(meta1,method="egger",k.min=1))# 调用函数并捕获输出
cat(output, file = "egger.txt", sep = "\n")

### Thompson####
metabias(meta1,method="Thompson",k.min=1)#线性图：发表偏倚
output <- capture.output(metabias(meta1,method="Thompson",k.min=1))# 调用函数并捕获输出
cat(output, file = "Thompson.txt", sep = "\n")

### linreg####
pdf(file="linreg.pdf", width = 8, height = 6)
metabias(meta1,method="linreg",plotit = TRUE,k.min=1)#线性图：发表偏倚
dev.off()
output <- capture.output(metabias(meta1,method="linreg",plotit = TRUE,k.min=1))# 调用函数并捕获输出
cat(output, file = "linreg.txt", sep = "\n")

